/** \file
 \brief Basic types and inline operations.
 */
#ifndef FASTMD_TYPES_INCLUDED
#define FASTMD_TYPES_INCLUDED 1

#include "yocto/bitwise.hpp"
#include "yocto/code/rand.hpp"

#include "yocto/memory/global.hpp"
#include "yocto/string.hpp"

#include <cmath>
#include <cstring>
#include <iostream>

using namespace yocto;

namespace FastMD 
{
	typedef memory::global         MemoryKind; //!< use new/delete
	typedef MemoryKind::allocator  Allocator;  //!< corresponding allocator
	
	//! Base class for objects..
	class Object : public object
	{
	public:
		explicit Object() throw() {}
		virtual ~Object() throw() {}
	private:
		YOCTO_DISABLE_COPY_AND_ASSIGN(Object);
	};
	
#define FASTMD_DOUBLE 0
	
#if FASTMD_DOUBLE == 0
	//==========================================================================
	// default: simple precision
	//==========================================================================
	typedef float Real;
#define FASTMD_F(X) X##f
#else
	//==========================================================================
	// FASTMD_DOUBLE: double precision
	//==========================================================================
	typedef double Real;
#define FASTMD_F(X) X
#endif
	
#define __sqrt      FASTMD_F(sqrt)
#define __cos       FASTMD_F(cos)
#define __sin       FASTMD_F(sin)
#define __pow       FASTMD_F(pow)
#define __floor     FASTMD_F(floor)
#define __ceil      FASTMD_F(ceil)
#define __log       FASTMD_F(log)
	
#define REAL_PI     FASTMD_F(3.14159265358979)
#define REAL_TWO_PI FASTMD_F(6.28318530717959)
#define __ONE       FASTMD_F(1.0)
#define __TWO       FASTMD_F(2.0)
#define __THREE     FASTMD_F(3.0)
#define __FOUR      FASTMD_F(4.0)
#define __HALF      FASTMD_F(0.5)
	
	//! ]0:1[ uniform random number
	inline Real   Rand() throw()  { return alea<Real>(); }
	
	//! [0:N-1] uniform random integer
	inline size_t LessThan(const size_t N) throw() { assert(N>0); return _rand.full<size_t>() % N; }	
	
	//! Normal(0,1) law
	Real Normal() throw();
	
	//! Initialize Random Numbers
	inline void   SRand() throw() { _rand.wseed(); }
	
	//! Fortran ANINT function
	inline Real   Anint( Real x ) throw() { return __floor( x + __HALF); }
	
	//! Periodic Boundary Condition
	/**
	 \param x  a position
	 \param L  the box length
	 \param iL 1/L
	 \return -L/2 <= x <= L/2
	 */
	inline Real   PBC1( const Real x, const Real L, const Real iL ) throw() { return x - L * Anint( iL * x); }
	
	//! vertex 3D
	class V3D
	{
	public:
		static const size_t SIZE;
		Real x;    //!< x coordinate.
		Real y;    //!< y coordinate.
		Real z;    //!< z coordinate.
		Real data; //!< mass, 1/mass,...
		
		//! default constructor.
		inline  V3D() throw() : x(0), y(0), z(0), data(0) {}
		
		//! construct from triplet.
		inline  V3D( const Real X, const Real Y, const Real Z ) throw() : x(X), y(Y), z(Z), data(0) {}
		
		//! default destructor.
		inline ~V3D() throw() {}
		
		//! copy constructor.
		inline  V3D( const V3D &other ) throw() : x( other.x ), y( other.y ), z( other.z ), data( other.data ) {}
		
		//! assign.
		inline  V3D & operator=( const V3D &other ) throw() { x=other.x; y=other.y; z=other.z; data=other.data; return *this;}
		
		//! vector addition
		inline void add( const V3D &rhs ) throw() { x += rhs.x; y += rhs.y; z += rhs.z; }
		
		//! vector subtraction
		inline void sub( const V3D &rhs ) throw() { x -= rhs.x; y -= rhs.y; z -= rhs.z; }
		
		//! vector multiplication.
		inline void mul( const Real a   ) throw() { x *= a; y *= a; z *= a; }
		
		//! this += a.rhs
		inline void muladd( const Real a, const V3D &rhs ) throw() { x += a*rhs.x; y += a*rhs.y; z += a*rhs.z; }
		
		//! this = a.rhs
		inline void mulset( const Real a, const V3D &rhs ) throw() { x  = a*rhs.x; y  = a*rhs.y; z  = a*rhs.z; }
		
		//! this -= a.rhs
		inline void mulsub( const Real a, const V3D &rhs ) throw() { x -= a*rhs.x; y -= a*rhs.y; z -= a*rhs.z; }
		
		//! dot product
		inline Real dot( const V3D &rhs ) const throw() { return (x * rhs.x) + (y * rhs.y) + (z * rhs.z); }
		
		//! norm, squared
		inline Real normsq() const throw() { return (x*x) + (y*y) + (z*z); }
		
		//! random vector with a given radius
		inline void rand( Real radius ) throw()
		{
			Real       u     = Rand();
			const Real theta = (REAL_TWO_PI) * Rand();
			u += u;
			u -= Real(1.0);
			const Real rho = radius * __sqrt( (Real(1.0) - u*u) );
			x =  rho * __cos( theta );
			y =  rho * __sin( theta );
			z =  radius * u;
		}
		
		//! current volume
		inline Real volume() const throw() { return x*y*z; }
		
		//! reciprocal vector
		inline void reciprocal() throw()
		{
			assert(x>0||x<0);
			assert(y>0||y<0);
			assert(z>0||z<0);
			x = __ONE/x;
			y = __ONE/y;
			z = __ONE/z;
		}
		
		//! special constructor: a.other.
		inline V3D( const Real a, const V3D &other ) throw() : x( a*other.x ), y( a*other.y ), z( a*other.z ), data(0) {}
		
		//! special constructor: b-a
		inline V3D( const V3D &a, const V3D &b ) throw() : x(b.x-a.x), y(b.y-a.y), z(b.z-a.z), data(0) {}
		
		//! set xyz to 0.
		inline void ldz() throw() { x=y=z=0; }
		
		//! I/O helper.
		inline friend std::ostream & operator<<( std::ostream &os, const V3D &v ) 
		{
			os << "[ " << v.x << " " << v.y << " " << v.z << " ]'";
			return os;
		}
	};
	
	//! Periodic Boundary Condition prototype
	typedef void (*PBC_Proc)( V3D &v, const V3D &box, const V3D &ibox);
	
	//! Fully Periodic
	inline void PBC_xyz( V3D &v, const V3D &box, const V3D &ibox) throw() {
		v.x = PBC1(v.x, box.x, ibox.x); 
		v.y = PBC1(v.y, box.y, ibox.y);
		v.z = PBC1(v.z, box.z, ibox.z);
	}
	
	//! No periodicity
	inline void PBC_none( V3D &, const V3D &, const V3D &) throw() {}
	
	//! units in 3D
	class U3D
	{
	public:
		size_t x,y,z;
		inline  U3D() throw() : x(0), y(0), z(0) {}
		inline  U3D( const size_t X, const size_t Y, const size_t Z ) throw() : x(X), y(Y), z(Z) {}
		inline ~U3D() throw() {}
		inline  U3D( const U3D &other ) throw() : x( other.x ), y( other.y ), z( other.z ) {}
		inline  U3D & operator=( const U3D &other ) throw() { x=other.x; y=other.y; z=other.z; return *this;}
		inline  size_t items() const throw() { return x*y*z; }
		inline  void ldz() throw() { x=y=z=0; }
		
		inline friend std::ostream & operator<<( std::ostream &os, const U3D &v )
		{
			os << "[ " << v.x << " " << v.y << " " << v.z << " ]'";
			return os;
		}
		
		inline friend bool operator==( const U3D &lhs, const U3D &rhs ) throw()
		{
			return (lhs.x==rhs.x) && (lhs.y==rhs.y) && (lhs.z==rhs.z);
		}
		
		inline friend bool operator!=(const U3D &lhs, const U3D &rhs ) throw()
		{
			return (lhs.x!=rhs.x) || (lhs.y!=rhs.y) || (lhs.z!=rhs.z);
		}
		
		inline size_t distance_to( const U3D &u ) const throw()
		{
			size_t ans = 0;
			ans += x >= u.x ? (x-u.x) : (u.x-x);
			ans += y >= u.y ? (y-u.y) : (u.y-y);
			ans += z >= u.z ? (z-u.z) : (u.z-z);
			return ans;
		}
	};
	
	
}

YOCTO_SUPPORT_C_STYLE_OPS(FastMD::V3D)
YOCTO_SUPPORT_C_STYLE_OPS(FastMD::U3D)

#endif
